﻿#include <iostream>
#include <fstream>
#include <math.h>
#include <exception>
#include <string>
#include <vector>
#include "func.h"
#include "parm.h"

//#include <stdlib.h>

using namespace std;




const int phi_bins = 100;


//输入文件名


int main(int argc, char** argv)
{
	
	double (*phi_t)[time_resolution] = new double[phi_bins][time_resolution]{};
	
	
	//自旋，内外半径，r2，环向宽度，输入输出文件名
	double a0 = stod(argv[1]);
	double rin = stod(argv[2]);
	double rout = stod(argv[3]);
	double r2 = stod(argv[4]);
	double dphcl = stod(argv[5]);
	ifstream input(argv[6]);
	ofstream output(argv[7]);

	dphcl = dphcl / 180 * Pi;

	double data[data_len] = {};

	double rcl0 = 0.5 * (rin + rout);
	const double omg = omega(a0, rcl0);
	const double period = 2 * Pi / omg;
	const double rclp = pow(rcl0, p0);

	try
	{
		while (true)
		{
			for (int i = 0; i < data_len; ++i) {
				if (!(input >> data[i])) throw runtime_error("read over");
			}
			//处理数据部分

			const double phi1 = phiem(data[xt], data[yt], data[sr], data[stheta]);
			//原团块对应的时间
			double time_origin_end = phi1 / omg;
			double time_origin_begin = time_origin_end - dphcl / omg;
			int time_origin_begin_int = normalization_time(time_origin_begin, period, time_resolution);
			int time_origin_end_int = normalization_time(time_origin_end, period, time_resolution);
			//像对应的时间
			double time_image_end = time_origin_end + data[str] + data[stth] - r2;
			double time_image_begin = time_image_end - dphcl / omg;
			int time_image_begin_int = normalization_time(time_image_begin, period, time_resolution);
			int time_image_end_int = normalization_time(time_image_end, period, time_resolution);
			//计算能段的强度

			int ji = round(elin * data[gt] / 0.1);
			double sg = data[ds] * pow(data[gt], 3) / rclp;


			//将符合时间要求(在第一个周期)的光子 加权加入结果中
			for (int t_image = time_image_begin_int; t_image < time_image_end_int; t_image++)
			{
				int t_in_period = t_image % time_resolution;
				if (t_in_period < 0)
				{
					t_in_period += time_resolution;
				}
				//一个是原团块的相位，一个是变化后团块的相位
				int phi_bin = floor((phi(data[xt], data[yt]) / Pi) * phi_bins/2);
				//int phi_bin = floor((phi1 /Pi) * phi_bins / 2) ;
				phi_bin = phi_bin % 100;
				if (phi_bin<0)
				{
					phi_bin += 100;
				}
				phi_t[phi_bin][t_in_period] += ji * sg;
			}
		}
	}

	catch (runtime_error err)
	{
		input.close();
	}



	for (int i = 0; i < phi_bins; i++)
	{
		for (int j = 0; j < time_resolution; j++)
		{
			output << phi_t[i][j] << " ";
		}
		output << endl;
	}
	output.close();

	return 0;
}
